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ABSTRACT 

A computational model aimed at predicting the flowfield of dump combustors is presented. 
The turbulent combustion model is based on the conserved scalar approach and on a 
convenient specification of its probability density function, which reduces the 
computation of the mean density to a closed form. Turbulence is modelled by means of the 
k- e model. The averaged conservation equations are solved by a technique based on a 
staggered grid and on the SIMPLE solver. 

The computational model is applied to a simple dump combustor to assess the computer time 
requirements and accuracy. The turbulent combustion model here introduced is shown to 
reduce the computer time by an order of magnitude when compared to evaluating the mean 
density by numerical quadrature. 


NOMENCLATURE 


a,b 

c j 

C 

C gl’ C g2 
C el> C c2 

f 

F ( f ) 
k 

h k (f ,f" 2 ) 

r\j 

H(f ,f" 2 ) 

J 

P 

_ A 

P 

P(f ) 


exponents of the beta -function 

coefficients of polynomial approximation of F(f) 
normalization constant of the beta function 

constants in the equation for the variance of the conserved scalar 

constants in the equation for the viscous dissipation rate 

constant in the expression for the turbulent viscosity 

conserved scalar 

quantity defined by eq. (4) 

turbulent kinetic energy 

factor defined by eq. (14) 

function defined by eq. (13) 
integral defined by eq. (11) 

order of the polynomial approximation of F(f) 
pressure 

quantity defined by eq. (2) 

probability density function of the conserved scalar 
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T 

u, v 
W 


x 

e 


$ 



a 


a f» a g ’ 


production terms defined by eqs. (3) 
radial coordinate 

universal gas constant 
temperature 

velocity components in x and r 
mixture mean molecular weight 
axial coordinate 
viscous dissipation rate 
generic quantity 
molecular viscosity 
turbulent viscosity 
density 

reference density 
molecular Prandtl number 

effective turbulent Prandtl numbers 


Subscripts 

( ) r derivative with respect to r 

( ) derivative with respect to x 


Superscripts 

(") conventional average 

f\j 

( ) Favre average 

( ) ' fluctuation about conventional average 

( ) " fluctuation about Favre average 


1 . INTRODUCTION 

Flameholding in subsonic combustion ramjets is achieved by using dump combustors, in 
which the flame is anchored to the chamber by a large region of recirculation. The 
ensuing elliptic character of the flow imposes however a considerable stress on the 
numerical simulation of such f lowfields . Special grid arrangements and computational 
techniques must therefore be used, and the computational requirements are demanding. The 
problem is compounded with that of providing an adequate representation of the effects of 
turbulent fluctuations on combustion. For nonpremixed flames the most popular approach 
is that making use of a presumed probability density function (pdf) of a conserved scalar 
quantity (ref. 1) . The effect of turbulent fluctuations is thus restricted to the mean 
density p, which is computed via an appropriate integral. The accurate evaluation of 
such integrals, which must be computed at each node point and at each computational step, 
requires a large number of intervals and takes therefore the largest share of the overall 
computing time. 

Three examples of previous applications of computer codes to the solution of dump 
combustor f lowfields in subsonic combustion ramjets are worthy of note. Harsha et al. 
(ref. 2) reduced the computer time to an acceptable level (4 to 15 minutes on a CDC 7600) 
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by separating the computation of the directed flow from that of the recirculation region: 
this requires an a priori knowledge of the reattachment length, so that this code cannot 
be considered to be a truly predictive tool. Further, no account was taken of the effect 
of fluctuations on combustion. Drummond (ref. 3) performed a unified calculation of the 
combustor flowfield, but the computer time was exceedingly high, even for nonreacting 
computations, in which 150000 steps and 4.6 hours on a CYBER 203 were needed. When 
reacting flowfields were considered the number of iterations rose to 200000, and the 
convergence of the solution was poor; the computer time is not reported. Also the 
combustion model did not account for the effects of fluctuations. The difficulty in 
obtaining a converged solution might be ascribed to the grid arrangement and the 
algebraic eddy viscosity model used, which appears not be suited to recirculating flows. 
A sounder computational model was presented by Vanka et al. (ref. 4) for ducted rockets, 
a ramjet variant. They correctly account for the effects of fluctuations by specifying 
the conserved scalar pdf as a battlement shaped function. The inherently 
three-dimensional character of the flowfield in this case requires using a rather coarse 
mesh (24 x 11 x 11) to limit the computing time to 25 minutes on an IBM 3033. 

The intent of the present paper is to show that, by introducing suitable assumptions to 
effect the thermochemical closure of the model, a very large reduction in computer time 
can be achieved. This is substantiated by applying the present model to a simple 
hydrogen-air dump combustor. In the following Section the describing equations are 
introduced, and then in Sec. 3 the proposed model for the thermochemical closure of such 
equations is presented. 


2. CONSERVATION EQUATIONS 

The bases of the analysis of the flowfield of nonpremixed turbulent combustion are the 
averaged conservation equations for mass, momentum and a scalar quantity f. The latter 
represents either a normalized enthalpy or atomic mass fraction: in fact, with the 

assumption that the kinetic energy is neglegible in the energy equation (implying 
relatively low- speed flow), that the diffusion can be represented in terms of a single 
coefficient, and that the Lewis number is unity, the equations for the above mentioned 
scalar quantities assume the same form, free of source terms. Favre- averaging (see e.g. 
ref. 5) is used to account for the effects of variable density: for a generic quantity 

$, it gives 

~ 

$ = — 

P 

where the tilde denotes the Favre mean, the overbar a conventional average, and the 
double prime the fluctuation about the Favre mean. Conventional averages are retained 
for the mean pressure p and density p . Turbulence is modeled according to the 

k - e model, which is generally accepted as the most reliable compromise between rather 
crude algebraic eddy viscosity models and second-order closure models. The latter are 
rather complicated and still not fully tested (see ref. 6, although considerable progress 
has been recently reported by Jones and Pascau, ref. 7) . To complete the model equations 
<v 

for the mean f and the variance f" of the conserved scalar (ref. 
1) are included in the formalism. Molecular diffusion terms are retained for 
completeness to account for near-wall and possible low Reynolds number effects (ref. 8), 
in particular at the inlet. As they are relatively important only in regions where the 
intensity of the fluctuations is low, such fluctuations are neglected in their 
evaluation. The steady, axisymmetric form of the conservation equations is considered. 
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Thus, the averaged equations for mass, axial and normal momentum, turbulent kinetic 

energy k , viscous dissipation rate e , and for the mean and the variance of the 
conserved scalar are 


P 1 f ~ 

puj =0 
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p is defined as 

P* “ P ‘ f {(P + P t ) [u x + y (r v) r ] + p k} 
and the production rates are given by 



[ rv 2 
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(le) 

(l f) 
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( 2 ) 


(3a) 

(3b) 


The turbulent viscosity is given as 


P- 


k 


r\j 

e 


a has the meaning of a molecular Prandtl number. The model constants are taken as 
(refs. 9 and 6) 


o ^ — 0.7 a k = 1 

C £i = 1.44 C e2 = 1.92 
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C, ( = 0.09 

Eqs . (1) represent a set of seven differential equations for the variables u, 

v, p*. f, k, 7 , f~ 2 . 

In order to close this set, the mean density p must be expressed as a function of the 
state variables considered. 


3. THERMOCHEMICAL CLOSURE 


Following the conserved scalar approach, all state quantities can be determined as a 
function of a single scalar quantity f, which is taken here as the elemental mass 
fraction of hydrogen. In fact, with the assumption that the pressure is thermochemically 
constant, the mixture composition and temperature (as well as any further state variable) 
are determined by chemical equilibrium computations , as a function of f alone . Thus from 
the equation of state: 


P 


P 


W 


T 


where R° is the universal gas constant and W is the mixture average molecular 
weight, the density too turns out to be a function of f alone. The above equation can be 
rewritten for convenience as 


P = 


ref 

F(f ) 


(4) 


where F(f) is a function of f determined by equilibrium computations; p is a 
reference density, introduced to make this function dimensionless. 


In order to evaluate the mean density, a probability density function (pdf) for the 
conserved scalar must be introduced, such that P(f) df is the fraction of time the 
conserved scalar is bounded between the values f and f+df. The mean density is then 
obtained from: 


P 


= P 


ref 


p(f) 
0 F(f) 


df 


(5) 


The pdf shape is assigned a priori, since it can be shown (ref. 1) that the result is 
relatively insensitive to the specific pdf assumed. Its parameters are determined on the 

basis of the mean f and the variance f" 2 of the conserved scalar, 
for which appropriately modelled differential equations are solved. In the present model 
(refs. 10, 11), developed after a suggestion by Libby (ref. 12), F(f) is approximated by 
a J-th order polynomial fit 


j 

F(f) 2 Y Cj f J 
j=o 


( 6 ) 


and a pdf is specified for the whole ratio P(f)/F(f) as a beta-function 


|||| = C f a ' 1 (l-f) b_1 (7) 

where C is a normalization constant and the exponents a and b of the beta -function are 
related to the mean and the intensity of the conserved scalar through the expressions 
(ref. 13) 
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r\j 2 r\j 

f (1-f) 

Aj 

f" 2 


Aj A*» r 

f (1-f)' 

Aj 

£||2 


- (1-f) 


( 8 ) 


Intermittency at the edges of the mixing layer is not accounted for explicitly in the 
present model. In any case, the beta- function itself accounts for intermittency to some 
extent, by giving a finite probability to the states f = 0 (if a < 1) and f = 1 (if 
b < 1). 


To determine the parameter C, the normalization condition 

i 

P(f) df = 1 
o 

is enforced. Its left-hand side can be written as 


1 (l-f) b _1 F(f)df 


Cj-tj (a, b ) 


j=o 


( 10 ) 


where use has been made of approximation (6); the integrals 1^ (j=0,l J) , defined 

as 


Ij(a,b) 


l 

f a " 1+J (l-f) b_1 df 

0 


(11) 


are expressed through the relation (derived after ref. 14) 


Ij (a,b) 


r (a+j ) r(b) 
r ( a+b+j ) 


Thanks to properties of the gamma- function, the summation appearing in (10) can be 
expressed as 
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j=o 
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This allows C to be determined. Upon substitution in (7) and then in (5), the mean 
density is obtained as 


P = 


ref 


H(f ,f" 2 ) 


with 
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where use has been made of (8) . 

This formulation has the advantage of reducing the integral appearing in (5) to a closed 
form. Further, the probability density function here introduced has been shown to give 
results fully consistent with experimental data (ref. 11). 

4. NUMERICAL ALGORITHM 

The present model for the mean density p has been introduced into the computer 
program JETFLO (ref. 15). This code is based on use of a staggered grid (ref. 16) and 
SIMPLE solver. The converged solution is obtained by iteration. Synthetic wall boundary 
conditions are enforced to avoid resolving the boundary layer, which would otherwise take 
the largest share of computing time. 

F(f) is approximated by a 12th order polynomial fit. The equilibrium conditions for the 
hydrogen-air flow considered here are computed by means of the equilibrium chemistry 
program STANJAN (ref. 17); the included species are Ar, H, H 2 , H 2 0, N, N 2 , NO, 
N0 2 , 0, 0 2 , OH. 


5. TEST CASE 

The present model is here applied to a simple cylindrical dump combustor. The main aim 
of this test case is to quantify the reduction of computing time achieved with the 
present turbulent combustion model, rather than to present an accurate computation of a 
practical combustor. Some limitations of the configuration being considered (i.e., the 
quite low Mach number and equivalence ratio) and of the computational implementations 
(i.e., the rather coarse mesh and the approximate specification of the inlet conditions) 
must be viewed in this perspective; they are introduced only for the sake of 
computational convenience and to limit the CPU time to an amount convenient for initial 
testing. The test case is computed with four different options, labelled 1 to 4: 

1) 28 x 28 grid, p given by the present model 

2) 15 x 15 grid, p given by the present model 

3) 15 x 15 grid, p computed by numerical quadrature over 200 intervals 

4) 15 x 15 grid, p computed by numerical quadrature over 100 intervals 

In runs 3 and 4 the probability density function P(f) of the conserved scalar is taken as 
a beta- function. 

The results of run 1 are presented in detail as the most accurate ones. A comparison of 
the results of runs 1 and 2 is used to give some indications as to what an extent the 
solution may be considered grid independent. The comparison of CPU time of runs 2 and 3 
is used to quantify the reduction achieved by the present model, as compared to 
evaluating p by numerical quadrature over 200 intervals. In order to assess whether 
200 intervals are adequate for estimating p , a comparison of the results of runs 3 
and 4 is used. 

The dump combustor considered in this test case is fed with hydrogen, injected through 
the central pipe (0.0077 m internal radius, 0.0154 m outer radius) at 1 m/s, and air, 
supplied through a 0.0115 m wide annulus at 10 m/s. The combustor internal radius is 
0.05 m, resulting in a rearward facing step 0.0231 m high at the inlet; its length is 0.3 
m. The velocity profiles at the inlet are assumed to be flat for each stream for 
simplicity. In any event this oversimplified assumption is unlikely to affect the overall 
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flowfield to a significant extent in elliptic flows and, moreover, the definition of the 
initial conditions for the test case being considered involves much larger incertainties, 
as will be apparent in the following paragraph. 

Figure 1 reports the main features of the combustor flowfield, as computed in run 1. In 
particular, Fig. la reports stream function isolines, i.e. , streamlines. A recirculation 
zone appears downstream of the dump plane and ensures flame stabilization. Due to the 
very low density of hydrogen, as well as to its low inlet velocity in this test case, 
such isolines are widely separated near the centerline, thus giving a rather poor 
representation of the flowfield. It is possible however to identify a smaller region of 
local recirculation near the hydrogen inlet, due to its inlet velocity being 
substantially lower than the other stream. It appears likely that such a recirculation 
spot reaches into the hydrogen supply pipe, thus invalidating the assumption of fully 
developed flow used at the inlet for computational purposes. This should however only 
affect the flowfield in the immediate neighbourhood of the hydrogen inlet, and have 
limited effects on the overall flow pattern. However this does suggest that a 
specification of the initial velocity profiles more accurate than that hypothesized above 
is unwarranted, at least in this case. The solution would be to extend the computational 
domain upstream into the hydrogen supply pipe. In Fig. lb the isolines of the logarithm 
of the mean value of the conserved scalar are shown. A rapid mixing is observed close to 
the plate separating the hydrogen and air streams. The stoichiometric combustion of 
hydrogen and air corresponds to f = 0.0284, so that the isoline -1.5 is very close to the 
mean location of the stoichiometric contour, albeit somewhat on the rich side. Of 
course, due to the fluctuations, such a front is highly oscillating. An estimate of the 
intensity of such oscillations can be obtained from Fig. lc which shows the isolines of 
the logarithm of the variance of the conserved scalar. The intensity of the fluctuations 
appears to be maximum close to the flame front. 

Figure 2 displays (for run 1) the isolines of the mean and the intensity of state 
quantities which are recovered from the mean and the variance of the conserved scalar via 
its assumed probability density function. For the mean density p, the form (12) is 
used. For any other state quantity (but the pressure) $, both conventional and Favre 
averages can be recovered as 


$ = 


l 

$(f) P(f) df 
o 


■'V 

$ 


p$ 

P 


1 

P 


l 

p(f) Hf) P(f) df 

o 


Further, the intensity of the fluctuations can be obtained as 


$ 


' 2 = 


l 


o 


[$(f)-4>] 2 P(f) df 


and similarly for $" 2 . 

In these expressions, $(f) denotes the equilibrium value of $ as a function of the 
conserved scalar. The evaluation of such integrals was achieved by numerical quadrature. 
However the contribution to the overall computing time is negligible as they need only be 
computed at the last iteration. The computation is further expedited by limiting the 
integration range to the domain 
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Fig. 1. Isoline plots: a) stream function (streamlines), b) logarithm of the mean of the 
conserved scalar, c) logarithm of the variance of the conserved scalar. 


631 










rather than [0,1]. 


Figures 2a through 2d show isolines of the average density, average temperature, rms of 
temperature fluctuations and average water vapour mass fraction, respectively. 
Conventional averages are reported for temperature as measurements by thermocouples give 
unweighted means, whereas Favre- averaging is used for the mass fraction as probe sampling 
results in close to density -weighted concentrations. It can be seen that the temperature 
at the exit of the combustor is rather low due to the large excess air involved in this 
particular test case. Mean temperatures about 2100 °K are predicted near the flame 
zone, with accompanying rms of fluctuations in excess of 770 °K. The mean water vapour 
mass fraction at the exit is rather low as well. 

We will now briefly compare the results of the different runs in terms of accuracy and 
CPU time . 

By comparing runs 1 and 2, it can be seen that the mean axial velocity component computed 
with the coaser grid (step size doubled in both directions) exhibits a maximum difference 
of 25% of the maximum velocity with respect to that computed with the finer grid. 
Although this may not represent a sufficiently grid- converged solution for most purposes, 
it is believed that the 28 X 28 grid gives results at least indicative of the overall 
flowfield configuration. Run 1 requires 443 iterations and 50.16 s CPU time on an IBM 
3090/600 processor; the computation is assumed to be converged when the residuals of all 

quantities are less than 10~ A . 

The comparison of the CPU time requirements of the present model ( p in closed form) 
and of a more conventional approach ( p evaluated by numerical quadrature) is 
performed on a 15 x 15 points grid, because the latter model involves very long 
computer runs. By comparing runs 2 (196 iterations and 6.11 s CPU time) and 3 (375 
iterations and 62.22 s) , it can be seen that the present model cuts the CPU time by a 
factor ten when compared to a numerical quadrature over 200 intervals . 

The question remain to be addressed whether 200 intervals are enough to ensure an 
accurate evaluation of p. To this ends, the results of run 3 can be compared to 
those of run 4 (numerical quadrature over 100 intervals) . The relative difference 
between the values of p computed with the two runs reaches a maximum of about 4%. 
This seems to indicate that the 200 intervals considered in run 3 might be adequate for a 
reasonably accurate computation, but 100 intervals might not. 

6. CONCLUSIONS 

The present model defines a soundly based treatment for nonpremixed turbulent combustion 
and, thanks to an appropriate assumption for the conserved scalar pdf, reduces the 
computer time requirements by about an order of magnitude in elliptic flows. 
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